green = "#009E73",
yellow = "#F0E442",
darkblue = "#0072B2",
red = "#D55E00",
purple = "#CC79A7")
# For nonsese...
#library(emo)
htmltools::tagList(rmarkdown::html_dependency_font_awesome())
library(tidyverse)
library(dplyr)
library(ggplot2)
theme_set(theme_bw())
load("IRdata.RData")
IRdata %>%
dplyr::count(vio, dem) %>%
IRdata %>%
dplyr::count(vio, dem)
IRdata %>%
dplyr::count(vio, dem) %>%
ggplot(mapping = aes(x =factor(vio), y = factor(dem))) +
geom_tile(mapping = aes(fill = n))+
labs(x = 'Violence', y = 'Regime type')
names(IRdata )
IRdata %>% dplyr::mutate(Z = pop*gdppc)
IRdata %>% dplyr::mutate(Z = pop*gdppc) %>% ggplot(aes(x =pop, y = gdppc)) + geom_tile(aes(fill = Z))
x <- LETTERS[1:20]
y <- paste0("var", seq(1,20))
data <- expand.grid(X=x, Y=y)
data$Z <- runif(400, 0, 5)
# Heatmap
ggplot(data, aes(X, Y, fill= Z)) +
geom_tile()
View(data)
IRdata %>% dplyr::mutate(Z = pop*gdppc) %>% ggplot(aes(x =pop, y = gdppc)) + geom_tile(aes(fill = Z)) +  scale_fill_continuous(guide = "colourbar")
library(correlation)
install.packages("correlation")
library(correlation)
data <- simulate_simpson(n=100, groups=10)
View(data)
library(ggplot2)
ggplot(data, aes(x=V1, y=V2)) +
geom_point() +
geom_smooth(colour="black", method="lm", se=FALSE) +
theme_classic()
#Simple correlation
correlation(data)
ggplot(data, aes(x=V1, y=V2)) +
geom_point(aes(colour=Group)) +
geom_smooth(aes(colour=Group), method="lm", se=FALSE) +
geom_smooth(colour="black", method="lm", se=FALSE) +
theme_classic()
correlation(data, multilevel = TRUE)
View(data)
simulate_simpson
##
library(dplyr)
data$Group
data2 <- data %>%
dplyr::mutate(team = ifelse(Group == "A" | Group == "B" | Group == "C", "Team A", NA)) %>%
dplyr::mutate(team = ifelse(Group == "D" | Group == "E" | Group == "F", "Team B", team))%>%
dplyr::mutate(team = ifelse(Group == "G" | Group == "H" | Group == "I" | Group == "J", "Team A", team))
View(data2)
table(data2$team)
data2 <- data %>%
dplyr::mutate(team = ifelse(Group == "A" | Group == "B" | Group == "C", "Team A", NA)) %>%
dplyr::mutate(team = ifelse(Group == "D" | Group == "E" | Group == "F", "Team B", team))%>%
dplyr::mutate(team = ifelse(Group == "G" | Group == "H" | Group == "I" | Group == "J", "Team C", team))
table(data2$team)
ggplot(data, aes(x=V1, y=V2)) +
geom_smooth(aes(colour=Group), method="lm", se=FALSE) +
facet_wrap(~team) +
theme_classic()
ggplot(data, aes(x=V1, y=V2)) +
geom_line(aes(colour=Group), method="lm", se=FALSE) +
facet_wrap(~team) +
theme_classic()
ggplot(data, aes(x=V1, y=V2)) +
geom_line(aes(colour=Group), method="lm", se=FALSE)
ggplot(data, aes(x=V1, y=V2, colour=Group)) +
geom_line()
ggplot(data, aes(x=V1, y=V2, colour=Group)) +
geom_line() +
facet_wrap(~ team)
ggplot(data2, aes(x=V1, y=V2, colour=Group)) +
geom_line() +
facet_wrap(~ team) +
theme_classic()
ggplot(data2, aes(x=V1, y=V2)) +
geom_smooth(aes(colour=Group), method="lm", se=FALSE) +
facet_wrap(~ team) +
theme_classic()
ggplot(data2, aes(x=V1, y=V2)) +
geom_smooth(aes(colour=Group), method="lm", se=FALSE) +
facet_wrap(~ team)
ggplot(data, aes(x=V1, y=V2)) +
geom_point(aes(colour=Group)) +
geom_smooth(aes(colour=Group), method="lm", se=FALSE) +
geom_smooth(colour="black", method="lm", se=FALSE) +
theme_classic()
ggplot(data, aes(x=V1, y=V2)) +
geom_point(aes(colour=Group)) +
geom_smooth(aes(colour=Group), method="lm", se=FALSE) +
geom_smooth(colour="black", method="lm", se=FALSE)
rm(list = ls())
library(bookdown)
install.packages("bookdown")
library(bookdown)
install.packages("modelsummary")
install.packages("likert")
help(likert)
??likert
library(likert)
help(likert)
data(pisaitems)
items29 <- pisaitems[,substr(names(pisaitems), 1,5) == 'ST25Q']
names(items29) <- c("Magazines", "Comic books", "Fiction",
"Non-fiction books", "Newspapers")
l29 <- likert(items29)
summary(l29)
plot(l29)
View(l29)
View(items29)
View(pisaitems)
View(items29)
summary(items29)
libs <- c('dplyr', 'tibble',      # wrangling
'stringr', 'readr',     # strings, input
'lubridate', 'tidyr',   # time, wrangling
'knitr', 'kableExtra',  # table styling
'ggplot2', 'viridis',   # visuals
'gganimate', 'sf',      # animations, maps
'ggthemes')             # visuals
invisible(lapply(libs, library, character.only = TRUE))
install.packages("kableExtra")
invisible(lapply(libs, library, character.only = TRUE))
## get data
infile <- "https://opendata.arcgis.com/datasets/dd4580c810204019a7b8eb3e0b329dd6_0.csv"
covid_de <- read_csv(infile, col_types = cols())
View(covid_de)
covid_de %>%
select(state = Bundesland,
county = Landkreis,
age_group = Altersgruppe,
gender = Geschlecht,
cases = AnzahlFall,
deaths = AnzahlTodesfall,
recovered = AnzahlGenesen,
date = Meldedatum)
covid_de %>%
select(state = Bundesland,
county = Landkreis,
age_group = Altersgruppe,
gender = Geschlecht,
cases = AnzahlFall,
deaths = AnzahlTodesfall,
recovered = AnzahlGenesen,
date = Meldedatum) %>%
mutate(date = date(date)) %>%
mutate(age_group = str_remove_all(age_group, "A"))
help(case_when)
covid_de %>%
select(state = Bundesland,
county = Landkreis,
age_group = Altersgruppe,
gender = Geschlecht,
cases = AnzahlFall,
deaths = AnzahlTodesfall,
recovered = AnzahlGenesen,
date = Meldedatum) %>%
mutate(date = date(date)) %>%
mutate(age_group = str_remove_all(age_group, "A")) %>%
mutate(age_group = case_when(
age_group == "unbekannt" ~ NA_character_,
age_group == "80+" ~ "80-99",
TRUE ~ age_group
)) %>%
mutate(gender = case_when(
gender == "W" ~ "F",
gender == "unbekannt" ~ NA_character_,
TRUE ~ gender
))
source("ISA2019/ggCirc.R")
library(postregplots)
library(tidyr)
library(ggplot2)
library(ggthemes)
library(viridis)
library(ggridges)
library(readr)
library(dplyr)
library(plyr)
library(data.table)
library(countrycode)
library("statnet")  # version 2016.9 with ergm 3.7.1
library("texreg")   # version 1.36.23
library("xergm")    # version 1.8.2 with btergm 1.9.0 and xergm.common 1.7.7
library(purrr)
library(amen)
library(abind)
library(plm)
load("Results/fit2_ongoing.RDta")
source("ISA2019/ggCirc.R")
specie <- c(rep("sorgho" , 3) , rep("poacee" , 3) , rep("banana" , 3) , rep("triticum" , 3) )
condition <- rep(c("normal" , "stress" , "Nitrogen") , 4)
value <- abs(rnorm(12 , 0 , 15))
data <- data.frame(specie,condition,value)
View(data)
ggplot(data, aes(y=value, x=specie)) +
geom_bar(position="stack", stat="identity")
library(ggplot2)
ggplot(data, aes(y=value, x=specie)) +
geom_bar(position="stack", stat="identity")
devtools::install_github("clauswilke/ggtextures")
install.packages("ggtextures")
help(barplot)
install.packages("randomcoloR")
3*40+4*49+3*47+10*41
rm(list = ls())
# set working directory
setwd("Replication_Data")
# Function to load packages
loadPkg=function(toLoad){
for(lib in toLoad){
if(! lib %in% installed.packages()[,1])
{install.packages(lib, repos='http://cran.rstudio.com/')}
suppressMessages( library(lib, character.only=TRUE))}}
# Load libraries
packs=c("maps", "rworldmap", "maptools", "sp", "spdep", "gstat", "MatchIt","cem",
"spatstat", "coefplot", "countrycode",  "cshapes", "arm", "stargazer","scales",
"texreg", "reshape2","ggplot2", 'foreign', 'car', 'lme4', 'dplyr', "ggmap","ggrepel",
"plotROC","pROC","purrr","tibble","magrittr","ROCR","caTools","xtable","rgdal",
"lattice","pscl","MASS")
loadPkg(packs)
Sys.info()
package_version()
remove(list = ls())
# required packages
library(ggplot2)
library(ggrepel)
library(zoo)
library(lme4)
library(dplyr)
library(scales)
library(ggpubr)
library(grid)
library(gridExtra)
install.packages(c("GADMTools", "tmap"))
library(ggplot2)
library(ggrepel)
library(zoo)
library(lme4)
library(dplyr)
library(scales)
library(ggpubr)
library(grid)
library(gridExtra)
# import Sciensano hospitalisations data
dat <- read.csv("https://epistat.sciensano.be/Data/COVID19BE_HOSP.csv")
View(dat)
# aggregate new intakes by province and date
dat <- aggregate(NEW_IN ~ DATE + PROVINCE, dat, sum)
View(dat)
# add new intakes for Belgium as a whole
belgium <- aggregate(NEW_IN ~ DATE, dat, sum)
View(belgium)
belgium$PROVINCE <- "Belgium"
col_order <- c("DATE", "PROVINCE", "NEW_IN")
belgium <- belgium[, col_order]
dat <- rbind(dat, belgium)
View(dat)
# transform date and provinces
dat$DATE <- as.Date(dat$DATE)
dat$PROVINCE <- factor(dat$PROVINCE,
levels = c(
"Antwerpen",
"BrabantWallon",
"Brussels",
"Hainaut",
"Li\xe8ge",
"Limburg",
"Luxembourg",
"Namur",
"OostVlaanderen",
"VlaamsBrabant",
"WestVlaanderen",
"Belgium"
),
labels = c(
"Antwerpen",
"Brabant Wallon",
"Brussels",
"Hainaut",
"Liège",
"Limburg",
"Luxembourg",
"Namur",
"Oost-Vlaanderen",
"Vlaams-Brabant",
"West-Vlaanderen",
"Belgique/België"
)
)
View(dat)
dat %>%
mutate(population = case_when(
PROVINCE == "Antwerpen" ~ 1857986,
PROVINCE == "Brabant Wallon" ~ 403599,
PROVINCE == "Brussels" ~ 1208542,
PROVINCE == "Hainaut" ~ 1344241,
PROVINCE == "Liège" ~ 1106992,
PROVINCE == "Limburg" ~ 874048,
PROVINCE == "Luxembourg" ~ 284638,
PROVINCE == "Namur" ~ 494325,
PROVINCE == "Oost-Vlaanderen" ~ 1515064,
PROVINCE == "Vlaams-Brabant" ~ 1146175,
PROVINCE == "West-Vlaanderen" ~ 1195796,
PROVINCE == "Belgique/België" ~ 11431406
))
View(belgium)
# compute NEW_IN by population size
dat <- dat %>%
mutate(population = case_when(
PROVINCE == "Antwerpen" ~ 1857986,
PROVINCE == "Brabant Wallon" ~ 403599,
PROVINCE == "Brussels" ~ 1208542,
PROVINCE == "Hainaut" ~ 1344241,
PROVINCE == "Liège" ~ 1106992,
PROVINCE == "Limburg" ~ 874048,
PROVINCE == "Luxembourg" ~ 284638,
PROVINCE == "Namur" ~ 494325,
PROVINCE == "Oost-Vlaanderen" ~ 1515064,
PROVINCE == "Vlaams-Brabant" ~ 1146175,
PROVINCE == "West-Vlaanderen" ~ 1195796,
PROVINCE == "Belgique/België" ~ 11431406
)) %>%
mutate(NEW_IN_divid = NEW_IN / population * 100000)
ggplot(
dat,
aes(x = DATE, y = NEW_IN_divid)
) +
geom_point(
size = 1L,
colour = "steelblue"
) +
labs(x = "", y = "Nombre d'hospitalisations (par 100,000 habitants) / Hospitalisaties (per 100,000 inwoners)") +
theme_minimal() +
facet_wrap(vars(PROVINCE),
scales = "free"
)
ggplot(
dat,
aes(x = DATE, y = NEW_IN_divid)
) +
geom_point(
size = 1L,
colour = "steelblue"
) +
labs(x = "", y = "Nombre d'hospitalisations (par 100,000 habitants) / Hospitalisaties (per 100,000 inwoners)") +
theme_minimal() +
facet_wrap(vars(PROVINCE),
scales = "free"
) +
geom_smooth(
se = FALSE,
col = "grey",
method = "gam",
formula = y ~ s(x)
)
ggplot(
dat,
aes(x = DATE, y = NEW_IN_divid)
) +
geom_point(
size = 1L,
colour = "steelblue"
) +
labs(x = "", y = "Nombre d'hospitalisations (par 100,000 habitants) / Hospitalisaties (per 100,000 inwoners)") +
theme_minimal() +
facet_wrap(vars(PROVINCE),
scales = "free"
) +
geom_smooth(
se = FALSE,
col = "grey",
method = "gam",
formula = y ~ s(x)
) +
geom_vline(
xintercept = as.Date("2020-05-04"), linetype = "dashed",
color = "darkgrey", size = 0.5
) +
geom_text(aes(x = as.Date("2020-05-04"), label = "1a", y = 9),
colour = "darkgrey", angle = 90, vjust = 1.3
)
# Create plot in dutch/fr
fig_trends <- ggplot(
dat,
aes(x = DATE, y = NEW_IN_divid)
) +
geom_point(
size = 1L,
colour = "steelblue"
) +
labs(x = "", y = "Nombre d'hospitalisations (par 100,000 habitants) / Hospitalisaties (per 100,000 inwoners)") +
theme_minimal() +
facet_wrap(vars(PROVINCE),
scales = "free"
) +
geom_smooth(
se = FALSE,
col = "grey",
method = "gam",
formula = y ~ s(x)
) +
geom_vline(
xintercept = as.Date("2020-05-04"), linetype = "dashed",
color = "darkgrey", size = 0.5
) +
geom_text(aes(x = as.Date("2020-05-04"), label = "1a", y = 9),
colour = "darkgrey", angle = 90, vjust = 1.3
) +
geom_vline(
xintercept = as.Date("2020-05-11"), linetype = "dashed",
color = "darkgrey", size = 0.5
) +
geom_text(aes(x = as.Date("2020-05-11"), label = "1b", y = 9),
colour = "darkgrey", angle = 90, vjust = 1.3
) +
geom_vline(
xintercept = as.Date("2020-05-18"), linetype = "dashed",
color = "darkgrey", size = 0.5
) +
geom_text(aes(x = as.Date("2020-05-18"), label = "2", y = 9),
colour = "darkgrey", angle = 90, vjust = 1.3
) +
annotate("rect",
ymin = -Inf, ymax = Inf,
xmin = as.Date("2020-03-15"), xmax = as.Date("2020-04-01"),
alpha = .2
) +
annotate("rect",
ymin = -Inf, ymax = Inf,
xmin = as.Date("2020-04-01"), xmax = as.Date("2020-05-01"),
alpha = .05
) +
annotate("rect",
ymin = -Inf, ymax = Inf,
xmin = as.Date("2020-05-01"), xmax = as.Date("2020-06-01"),
alpha = .2
) +
annotate("rect",
ymin = -Inf, ymax = Inf,
xmin = as.Date("2020-06-01"), xmax = as.Date("2020-06-05"),
alpha = .05
) +
labs(
title = "Evolution des admissions hospitalières / Evolutie van de hospitalisaties - COVID-19"
) +
scale_y_continuous(breaks = seq(from = 0, to = 10, by = 2), limits = c(0, 10)) +
scale_x_date(labels = date_format("%m-%Y"))
fig_trends
## adjust caption at the end of the trend figure
caption <- grobTree(
textGrob(" * Lignes solides : courbes ajustées aux observations / Volle lijnen : gefitte curves \n * Lignes pointillées : phases de déconfinement 1a, 1b & 2 / Gestippelde lijnen: fases afbouw lockdown maatregelen 1a, 1b & 2",
x = 0, hjust = 0, vjust = 0,
gp = gpar(col = "darkgray", fontsize = 7, lineheight = 1.2)
),
textGrob("Niko Speybroeck (@NikoSpeybroeck), Antoine Soetewey (@statsandr) & Angel Rosas (@arosas_aguirre) \n Data: https://epistat.wiv-isp.be/covid/  ",
x = 1, hjust = 1, vjust = 0,
gp = gpar(col = "black", fontsize = 7.5, lineheight = 1.2)
),
cl = "ann"
)
caption
library(GADMTools)
library(RColorBrewer)
library(tmap)
## sf structure
map <- gadm_sf_loadCountries(c("BEL"), level = 2, basefile = "./")$sf
rm(list = ls())
rm(list=ls())
library(lme4)
library(arm)
library(plm)
library(dplyr)
library(texreg)
library(extrafont)
library(ggplot2)
library(showtext)
library(ggthemes)
library(scales)
font_add('KaiTi', 'KaiTi.ttf')
showtext_auto(enable = TRUE)
######### fixed-effect model with SD clustered by ccode
setwd("~/Documents/Research/SwapProject/PC_Paper/PC_Code/Replication")
source("RFunctions_modeling.R")
# Note: all Dvs have been led for one period (which is equivalent to lagging all IVs )
load("data.RData")
## Let's call the models mixed-effect models (economists like it more)
summary(data$year)
data <- data %>%
dplyr::filter(ccode !=710 )
data <- data %>%
dplyr::mutate(defense_usa_lag1 = ifelse(is.na(defense_usa_lag1), 0, defense_usa_lag1)) %>%
dplyr::filter(year < 2019 & year > 2008)
# m0: simple model
m0 <- MLMs_fit(dvs = c('distance_idealpoint_chnlog'),
ivs = c('distance_idealpoint_chnlog_lag1','swap_dummy_lag1','distance_idealpoint_usalog_lag1','minidist_usalog_lag1',
'defense_usa_lag1'),
modname = c('base'),
data = data)
summary(m0$base)
